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SUMMARY 


A three-node flat shell element with C° rotation fields has been developed for analysis of arbitrary composite 
shells. The element may consist of any number of orthotropic layers, each layer having different material properties 
and angular orientation. The formulation includes coupling between bending and extension, which is essential for 
analysis of unsymmetric laminates. Shearing deflections are included, since laminated and sandwich construction 
frequently results in shear stiffness much smaller than bending stiffness. Formulation of the element is straight- 
forward, and calculation of its stiffness matrix is simple and fast. Convergence of solutions with mesh refinement 
is uniform for both thin and thick shells and is insensitive to element shape, although not as rapid as some other 
elements that lack one or more capabilities of the newly developed element. An experimental verification of the 
shell element is reported in the appendix. 

INTRODUCTION 

Significant improvements in aircraft flutter speeds and other performance characteristics are possible through cre- 
ative use of composite materials, as has been pointed out by Weisshaar (1987). Finite-element analysis of structures 
of this type requires shell elements representing layered composites that may be unsymmetric and may frequently 
use honeycomb or foam cores. The element described here was developed for this purpose. 

Since the beginning of research on the finite-element method, a great deal of effort has been expended on devel- 
opment of shell elements. Many papers on this subject have been published (Gallagher, 1975; 1978). However, the 
perfect shell element has yet to be invented. 

Triangular flat-elements having displacements and rotations at the comer nodes as degrees of freedom (dof) are 
appealing for practical reasons. They can easily model arbitrary shell geometries with general supports, cutouts, and 
beam stiffeners. These elements have a total of 18 dof (three translations and three rotations at each node) or 15 dof 
(three translations and two rotations), depending on whether the rotation about the normal, which has zero stiffness, 
is included as a degree-of-freedom. 

One of the earliest shell elements was a flat triangular element developed by Melosh (1966) and Martin (1967). 
Versions of this element were used in several early programs including ELAS, SAMIS, and DYNAL. At that time, 
it was the only element that used the transverse shear modulus, and, consequently, it was the only element suitable 
for sandwich plates and shells. In 1971, Utku (1973) developed a new version with a straightforward treatment of 
shearing stiffness which is theoretically much improved. 

The element developed here follows the derivation by Utku (1973). However, integration through the thickness 
is done stepwise, required for representing laminates, and two additional terms of the energy integrals are retained, 
since they do not vanish for unsymmetric laminates. 

This element converges at about the same rate for thick or thin plates, regardless of element shape. It appears to be 
quite robust and is capable of correctly modeling behavior of general, unsymmetric, sandwich, laminated composite 
shells. Convergence is slower than that of some other three-node elements that do not include shearing deflections. 
Because of its straightforward formulation, this element may provide a basis for development of a dynamic element 
(Gupta, 1979), which could result in a significant improvement of numerical efficiency in analysis of problems of 
aircraft vibration and flutter. 

ELEMENT FORMULATION 

The triangular element is developed in a local coordinate system with origin at the middle surface centroid, and 
the x axis is parallel to the element 1-2 edge, as shown in figure 1. The z axis is perpendicular to the plane of the 



element. In the computer implementation of the element, material properties for a layer are defined with respect to 
a material coordinate axis which may be oriented at an angle a measured from the element x axis to the material 
I axis. The material 3 axis coincides with the element 2 axis. Displacements in the x, y, z directions are denoted 
by u,v,w. The usual right-hand vector convention is used for directions of rotations, which are denoted by 6 X , 
6 y , and d 2 . The rotations 0 X and 0 V are rotations of lines originally perpendicular to the middle surface of the 
undeformed element. 

Using the comma convention for differentiation, the midsurface strains e 0 are given by 


and curvatures are given by 


e a = 


w o,x 

v o t y 

u o,y + v o,x 



— 9y tX 


* ^!ix 

X = 

^x,y 

= 

w ',yy 


$x,x “ 8y,y 


. 2 w \*v . 


In this derivation, the total deflection w is in the z direction 


U) = yj f + yj* 

where w f is the Kirchhoffian part of the deflection and w* is the part due to out-of-plane shear. 
At any point in the element, the strains are given by 

e = e 0 — zX 


and stresses are given by 

a = De = De 0 — DzX 

where I? is a 3 by 3 material property matrix for plane stress in the lamina. 
The out-of-plane shear strains are given by 
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and the corresponding shear stresses are given by 
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The shear moduli in D' are reduced by a factor of 5/6 to account for nonuniform distribution of shear stress, as is 
done in many other programs. 

In this formulation, a straight line, which is normal to the middle surface before loading, remains straight, but 
not normal to the deformed middle surface after loading. The slopes of the middle surface are w jX and w >yi whereas 
the slopes of the line that was originally normal to the middle surface become ~9 y = w iX — and 9 X = w iV — w* 
and 9 X = — v iZ and 9 y - u tZ . 

The total strain energy in the element is given by 

U - — J € T adAdz + - - J rdAdz 
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The energy expression can be expanded into five terms, as follows: 

U= ^ J (?De 0 dAdz+ j j X T DXz 2 dAdz + j J D'^dAdz + j J X T De 0 zdAdz+ j J <?DXzdAdz 

The last two integrals represent coupling between middle-surface deformations and bending or twisting which vanish 
if the laminate is symmetric with respect to the middle surface, but must be retained in this development. 

In integrating through the thickness, it is assumed that each layer is a different homogeneous orthotropic material, 
but none of the matrices depends on 2 . For example, the last integral in the aforementioned energy expression 
consequently becomes 

[ £ DXdA 

4 .=1 Ja < 

where hi is the 2 coordinate of the bottom of layer i and n is the number of layers. 

The element matrix is obtained by differentiation of the energy integral with respect to the nodal displacements, 
as is usually done 

K = J^L 

11 duiduj 

where U{, Uj represent all of the nodal displacements and rotations. 

To proceed with development of the element, it is necessary to define the matrices that relate strains in the element 
to nodal displacements. These expressions are taken from the paper by Utku (1973), and the details of the derivation 
will not be repeated here. The nodal displacements in the x direction at the three nodal points are represented by 


u = 


u i 
U2 
U 3 


and the same convention is used to represent the other components of nodal displacement. 

The displacements in the interior of the element u , v , 6 X , 8 y , and w* are each interpolated from their nodal values 
by the same linear function. For example, 


u = 


2 A 


[x y 1] 
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X2 V3 — X3 V2 


t/3 - y \ 
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[u] 


To get the interpolation for deflection in the y direction, u is replaced by v and u by v, and so forth. In this 
expression, A represents area of the element, and x\ and y\ are the coordinates of node 1, and so forth. 

It follows that the strains on the middle surface e Q are given by 


2 a [MN] 
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and the curvatures X are given by 


where 
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and 


M = 


yi - 1/3 yj - i/i 1/1 - 1/2 
0 0 0 
X3 — X2 Xl — X'i X2 — Xl 


Development of the expression for out-of-plane shear strain 7 in terms of nodal displacements is the principal 
contribution of Utku’s (1973) paper. It is given by 
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Perhaps it should be noted that some terms are missing from the matrices which are the end result in Utku (1973). 

In the computer implementation of this element, the matrices relating strains to nodal displacements are calcu- 
lated from the aforementioned algebraic expressions. The remaining problems are calculated numerically, layer by 
layer. Since the strain-displacement matrices contain only constants, integration over the area is trivial. Each of 
the five energy integrals is processed to yield the corresponding submatrix of the element, stiffness matrix, and the 
results are added into the proper location in the element stiffness matrix. The program then goes on to process the 
next layer which may have different material properties. The first term of the energy integral, for example, represents 
the membrane deformation part of the strain energy U m . 

U m = \ J ^Ue 0 dAdz 

Integrating through the thickness, this becomes 


1 n r 

U m = 1 - hd / ? 0 De 0 dA 

1 i=i Ja > 


where hi is the z coordinate of the bottom of layer i and ra is the number of layers. 

Substituting for strains in terms of nodal displacements and taking second partial derivatives of the energy with 
respect to nodal displacements, results in the submatrix k m , which represents the membrane part of the element 
stiffness matrix. 

*” - ~ WMNfDdM N] 

Using the same process 


k>b 

k a = 


1 

12A 

1 

16 A 


- h blN -MfDdN- M] 

i= 1 
n 

-hdlHfD'dH] 
i= 1 


4 


kAl = A £( hi 1 - ) [ M NfD[ N-M ] 

i=i 

= nD^i - -A/] r A[Af iV] 

,=1 

Rotations 0 X and 8 y are row degrees of freedom for k^\ • Displacements u amd t> are row degrees of freedom for 
kAi . The element stiffness matrix K is 


K — k m + k), + k s + kAi + kA2 

One term in each of the aforementioned summations, which represents stiffness contribution of the current layer, 
is computed. Each stiffness coefficient from each submatrix is added into the appropriate location in the element 
stiffness matrix. The program then goes on to the next layer until all layers have been processed. The element 
stiffness matrix is then complete. 

ELEMENT PERFORMANCE 

The three-node shell element developed herein should prove convenient for modeling irregular shapes with 
cutouts, as well as for analysis of unsymmetric sandwich composite laminates and includes out-of-plane shear de- 
formation. In comparison with other elements, it may be termed LU71 (Laminated-Utku-1971). Calculations show 
that the LU71 gives accurate solutions for both thick and thin plates, converges uniformly for various shapes, and 
appears robust. Convergence with mesh refinement is slower than that of the two elements used for comparison, but 
neither of those elements can fill the need for which the LU71 was developed. 

Figure 2 shows one of the test problems used to verify the element, which illustrates the coupling of bending 
and in-plane deflection. Table 1 shows calculated deflections for four cases, all with the same loading. These results 
are exact and can be easily calculated by hand because strains are constant throughout. The first case is a four-layer 
symmetric laminate with fibers only in the 0° and 90° directions. The second case represents modeling where 0° and 
90° elements are overlayed. Results are identical. The third case is an unsymmetric laminate. Loads applied at the 
middle surface cause bending and out-of-plane deflection. The fourth case is identical to the third, except that the 
laminate thickness is 0.1 in. instead of the 1.0 in. used in the first three cases. Deflection at the center of stiffness 
is the same in all four cases. The increase in the value of the deformation in the X direction UX for the third and 
fourth cases is due to rotation and the calculation of the UX at the middle surface. 

In this example, Material I was used as defined by Reddy (1980) with E\ = 10 .6 x 10 6 psi. 

E\ I Ei = 25 
Gn/Ei = 0.5 
Gii / E 2 = 0 .2 
v \2 = 0 .25 


Figure 3 provides a convergence comparison of the LU71 element with the DKT element (Batoz and others, 
1980), which is type 53 in ANSYS (DeSalvo and Swanson, 1985), and with Reddy’s eight-node isoparametric 
element (Reddy, 1980) for a thin simply supported plate. All converge to the same result. The DKT type B solution 
has element diagonals at —45° rather than +45°, and convergence is slower, but still very good. 

The convergence in figure 4 is compared with that in figure 3 for a moderately thick plate. The DKT element 
is much too stiff because it has no shearing deflections. The LU71 element converges to the same result as that of 
Reddy’s (1980) isoparametric element. 
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In figures 3 and 4, Reddy’s result is plotted as being equivalent to 32 triangular shell elements, since it is estimated 
that 32 triangular shell elements represent about the same computational effort as the four isoparametric quadrilateral 
elements actually used. Reddy uses 2 by 2 integration in his examples, but for some 8-node elements, this can 
result in an ill-conditioned or even singular matrix in structures with no or few support constaints. In application to 
airframes, it is believed that 3 by 3 integration would be required. With a 3 by 3 integration, the result will be a little 
less accurate than that shown in the figures. 

The DKT element (Batoz and others, 1980) gives outstanding performance in most bending problems. Its con- 
vergence rate is, however, sensitive to boundary conditions and element shape. For membrane deformations, it is 
a constant stress triangle, as in the LU71. Lack of shearing deformations is its greatest weakness in analysis of 
composites. 

Utku (1973) gives seven examples comparing the performance of the shell element formulation used here with 
an earlier formulation for isotropic homogeneous plates. The examples show uniform convergence regardless of 
element shape. The formulation used here converges faster than the earlier one for equilateral triangular elements; 
but the convergence is slower when the element has an obtuse angle. 

Table 2 shows some results for thin plates obtained by using the LU71. The example is a 48-in. 2 homogeneous 
isotropic plate with E = 10.6 x 10 6 psi, v = 1/3 and a concentrated load of 16 lb at the center. All edges were 
fixed. One quarter of the plate was modeled by 242 elements in a regular mesh. Elements were right triangles with 
the short side 2.18 in. long. Since the LU71 includes shear deformation, it should give larger deflections for the 
thicker plates but the same deflection for very thin plates. Table 2 shows that the LU71 gives results as expected at 
the plate side length/thickness ratio a/h of 500. At a/h = 1000 , the LU71 has a small error, for thinner plates, the 
errors become greater. 

CONCLUDING REMARKS 

The three-node shell element proves to be convenient for modeling irregular shapes with cutouts, and is also 
suitable for analysis of unsymmetric sandwich composite laminates. It appears to be robust and gives accurate 
results for both thick and thin shells. 

Formulation of the element is straightforward. It is believed that this formulation can be used as a basis for a 
dynamic element that may provide much improved numerical efficiency. 


Ames Research Center 
Dryden Flight Research Facility 
National Aeronautics and Space Administration 
Edwards, California, March 1, 1988 
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APPENDIX— COMPARISON OF EXPERIMENTAL DEFLECTIONS OF 
GRAPHITE COMPOSITE PLATES WITH FINITE ELEMENT CALCULATIONS 


Introduction 

Three graphite-composite plates were built and tested to verify the performance of the LU71 finite element. The 
LU7 1 element was developed during 1986 and implemented for testing in a program called LCSHEL during 1987. It 
is a three-node C° element for analysis of plates and shells made of composite materials with an arbitrary symmetric 
or unsymmetric stacking sequence. The LU71 was developed to meet an immediate need in the STARS program 
(Gupta, 1984) for analysis of airframes made of unsymmetric laminated composites and also to provide a foundation 
for future development of a finite dynamic element with greatly improved numerical efficiency. The experimental 
work described here was intended to verify the treatment and manipulation of material properties for symmetric and 
unsymmetric laminates and the overall correctness of analysis by using the element. 

Plate Fabrication 

Three graphite fiber composite plates were fabricated. The graphite fiber was supplied in a 12-in.-wide roll 
of unidirectional fiber with a few cross-strands of glass fiber to maintain shape and fiber spacing. The specified 
weight of this material was 0.033 lb/ft 2 , and specified tensile strength of the fiber was 450,000 psi. The resin used 
was Hexcell 2410 with 2184 hardener. Plates were placed on a flat surface, vacuum bagged, and cured at room 
temperature. 

Plate A was a four-layer unidirectional plate 9.25 by 12.0 in. with an average thickness of 0.0303 in. The fiber 
fraction was 0.572 by weight. Calculated volume fractions were Vf = 0.474 fiber, v m = 0.458 resin, and 0.068 
voids. This plate was cut into strips for tension tests in the fiber and transverse directions for determination of 
material properties. 

Plate B was a six-layer unsymmetric (0/0/0/90/90/90) plate 12.0 by 12.0 in. with an average thickness of 
0.04311 in. The fiber fraction was 0.582 by weight. Calculated volume fractions were v/ = 0.494 fiber, v m = 0.480 
resin, and 0.0259 voids. Thickness was measured at 1-in. grid points on a surface plate with a 1/10,000-in. dial in- 
dicator. The average thickness was 0.04311 in. with a standard deviation of 0.0058 in. This plate was trimmed to 
11.3 by 11.3 in. for use in a five-point bend test. 

Plate C was an eight-layer symmetric (0/0/90/90/90/90/0/0) plate, 12.0 by 12.0 in. with an average thickness of 
0.0549 in. The fiber fraction was 0.637 by weight. Calculated volume fractions were Vf = 0 .517 fiber, v m = 0 .399 
resin, and 0.084 voids. Thickness was measured at 1-in. grid points on a surface plate with a 1/10,000-in. dial 
indicator. The average thickness was 0.0549 in. with a standard deviation of 0.0046 in. This plate was trimmed to 
1 1.3 by 11.3 in. for use in a five-point test. 

Material Properties 

Material properties were determined from tests on specimens cut from plate A. Specimens were approximately 
1 in. wide by 9 in. long. Metal plates 2.5 in. long were glued to the ends of the specimens (to avoid crushing), leaving 
a 4-in. test length. Biaxial electric resistance strain gages were glued to both sides of the specimens and all gages 
were read individually. 

Material properties in the fiber direction were consistent. Unfortunately, the test machine for measuring trans- 
verse properties failed to work properly, and some data were lost. Since there was some uncertainty in the correct 
value of the transverse modulus, the material properties in table 3 show a high estimate for E 2 (upper bound) and a 
low estimate for Ei (lower bound). 
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Finite-element calculations were performed for both sets of properties so that the effect of this uncertainty could 
be assessed. Adjustments in properties to account for differences in fiber volume fractions were done using relations 
from Chamis (1984) and Jones (1975). 

The in-plane shear modulus Gn was determined from a four-point bend test, supported at three comers and 
loaded at the fourth comer. The shear modulus is determined directly as given by Chandra (1976). When limited to 
the linear range and neglecting self- weight, this simplifies to 


G n 


3 q 2 P 
h 2 w 


in which a is plate side length, h is plate thickness, P is load, and w is deflection at the point of loading. The 
value of Gn is larger than would be expected from Chamis (1984), but numerical experiments show that this has a 
negligible influence on results. This larger value of in-plane shear modulus could be associated with slight waviness 
in the fibers. 


Figure 5 shows the subscripting convention and notation for elastic constants. 


Plate Bending Tests 

Plates B and C were tested in five-point bending with supports at each of the four comers and a concentrated load 
at the center (fig. 6). Supports were set 0.35 in. from the comers of the plates so that the distance between supports 
was 10.6 in. Loading was by deadweights, and deflections were measured by a microscope with a resolution of 
1/10,000 in. 

Plots of load versus deflection data are shown for the six-layer unsymmetric plate in figure 7 and for the eight- 
layer symmetric plate in figure 8. Lines through the data points were fitted by least squares. The plots show good 
linearity for both loading and unloading. 

Comparisons to Calculations 

The finite-element grid used for calculations is shown in figure 9. Calculations were done with the LU71 ele- 
ment in the LCSHEL program. Additional calculations, not included here, were done with ANSYS, and agreement 
between the two programs was good. 

Experimental and calculated results are compared in table 4. Since there was some uncertainty about the values 
of ^.calculations were done separately using the high estimate (upper bound) for E 2 and the low estimate (lower 
bound) for E 2 . 

Results for the six-layer unsymmetric plate (plate B) show excellent agreement. Using the upper-bound prop- 
erties, the calculated result is 2.3 percent stiffer than the experimental result; using the lower-bound properties, the 
calculated result is 2.2 percent softer than the experimental result. 

Agreement between results from experiment and calculations for the eight-layer symmetric plate is not optimum, 
but still adequate. Using the upper-bound properties results in a calculation 5.7 percent stiffer than the experiment, 
while using the lower-bound properties results in a calculation that is 2.2 percent stiffer than the experiment. The 
average deviation between calculations and experiment is 3.9 percent. 


Conclusion 


Comparison of results from plate bending experiments and finite-element calculations show very good agree- 
ment, and support the conclusion that the finite-element program is using composite material properties correctly. 
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It is believed that the resulting discrepancies are due primarily to imperfections in technique for fabricating the 

composite plates and performing basic property measurements. 
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TABLE 1. DEFLECTIONS FOR IN-PLANE LOADING 


Stack 

sequence 

ux, 

10- 4 in. 

UY, 

10 ~ 6 in. 

UZ, 

10“ 2 in. 

10 ~ 3 rad 

0/90/90/0 

-0.36207 

0.69629 

0 

0 

Superimposed 
0 and 90 

-0.36207 

0.69629 

0 

0 

0/90/90/0 

-1.0038 

1 .9303 

0.13898 

-0.27797 

0/0/90/90 

-10.038 

19.303 

13.898 

-27.797 


TABLE 2. CENTER DEFLECTION OF THIN SQUARE PLATES WITH FIXED EDGES 3 


h, in. 

a/h 

b/h 

Thin plate 
theory, in. 

LU71, 

in. 

Difference, 

percent 

0.96 

50 

2.27 

0.2180(10~ J ) 

0.2329 x 10 

+6.8 

0.48 

100 

4.55 

0.1744(10- 2 ) 

0.1840 x 10 —2 

+5.5 

0.096 

500 

22.7 

0.2180 

0.2231 

+2.3 

0.048 

1000 

45.5 

1.744 

1.684 

-3.4 

0.032 

1500 

68.2 

5.887 

5.287 

-10.0 

0.024 

2000 

90.9 

13.95 

11.56 

- 17.2 


“Plate side length = o = 48 in. 

Thickness = h 

Length of short element side = b = 2 .18 in. 


TABLE 3. MATERIAL PROPERTIES 



Low Ea 
estimate 

High Ei 
estimate 

Fiber modulus, Ef, 10 6 psi 

30.13 

30.09 

Resin modulus, E mt 10 6 psi 

0.2138 

0.2484 

Vm 

0.3000 

0.3000 

v s 

0.1383 

0.0811 

Plate B, six layers 

Ei, 10 6 psi 

15.38 

15.38 

£ 2 , 10 6 psi 

0.8463 

0.9784 

^12 

0.2181 

0.1886 

^21 

0.0120 

0.0120 

G 12 , 10 6 psi 

1.849 

1.849 

Plate C, eight layers 

E\ , 10 6 psi 

15.64 

15.65 

Ei, 10 6 psi 

0.7231 

0.8365 

^12 

0.1911 

0.1868 


0.0088 

0.0100 

Gn, 10 6 psi 

1.938 

1.938 



TABLE 4. COMPARISON OF EXPERIMENTAL MEASUREMENT AND FINITE-ELEMENT 
CALCULATIONS OF STIFFNESS (LOAD/DEFLECTION RATIO) OF GRAPHITE COMPOSITE PLATES 


Plate 

Properties 

Finite-element 
calculated stiffness, 
lb/in. 

Experimental 

stiffness, 

lb/in. 

Difference, 

percent 

Plate B 

(Six-layer) 

High Ei estimate 

6.504 

6.359 

+2.3 


Low Ei estimate 

6.223 

6.359 

-2 .2 


Average 

6.364 

6.359 

0.0 

Plate C 

(Eight-layer) 

High Ei estimate 

14.45 

13.67 

+5.7 


Low Ei estimate 

13.97 

13.67 

+2 .2 


Average 

14.21 

13.67 

+3.9 



Figure 1. Coordinate systems. 


t UY 1 0Y 
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X 

Material I 




Figure 2. In-plane loading example. 
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Stack sequence 
[0/90/90/0] O 

a = 20 in. 

h = 0.2 in. 

a/h = 100 


Reddy’s 2 by 2 element 
LU71 program 

ANSYS DKT element type B 
ANSYS DKT element type A 



Stack sequence 
[0/90/90/0] 

a = 20 in. O Reddy’s 2 by 2 element 

h = 1 in. LU71 program 

a/h = 20 ANSYS DKT element 



Figure 3. Comparison of convergence for thin plate. Figure 4. Comparison of convergence for thick plate. 


The strains for plane stress due to 
load in direction 1 are 



where the loading direction is 
denoted by left superscript. 
Similarly, the strains due to load 
in direction 2 are 
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therefore, the stiffness matrix for 
the stress-strain relations for plane 
stress in orthotropic material in 
terms of the engineering constants 
are 
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Figure 5. Material property definitions (Jones, 1975, pp. 39-46). 
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Figure 7. Bending test of six-layer unsymmetric plate. Figure 8. Bending test of eight-layer symmetric plate. 



Figure 9. Finite-element mesh used for calculation of stiffness of graphite composite plates. 
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